Rocky Shore Tanaids Compostion

library(readxl)
library(reshape2)
library(doBy)
library(dplyr)
library(stringr)
library(ggplot2)
library(ggforce)
library(vegan)
library(ggdendro)
library(patchwork)
library(ggrepel)
library(rjson)
library(car)
library(plyr)
library(pheatmap)
library(RColorBrewer)

ggplot theme setting

large <- theme(legend.title = element_text(size=15),
        legend.text = element_text(size=15),
        axis.title = element_text(size=15),
        axis.text = element_text(size=15),
        strip.text = element_text(size=15))

rotate <- theme(axis.text.x = element_text(angle = 90, hjust=1))

no_strip <- theme(strip.background = element_rect(colour=NA, fill=NA),
                  strip.text = element_text(colour=NA))

Preparing tanaid data

ta0 <- read_excel("../Data/Tanaidacea community.xlsx", sheet = 3)
ta <- ta0[, 1:2]
# Create factor "Microhabitat"
ta$Microhabitat <- str_replace(ta0$Microhabitat, " \\s*\\([^\\)]+\\)", "")
ta$Microhabitat2 <- ta$Microhabitat
ta$Microhabitat2[!(ta$Microhabitat2=="Gravel, sand or silt" | ta$Microhabitat2=="Tube of Eunice taoi")] <- "Macroalgae"
# Create factor "Replicate"
ta$Replicate <- str_extract(ta0$Microhabitat, "(?<=\\()\\d+(?=\\))")
ta <- cbind(ta, ta0[,-1:-3])
# Remove zero rows
ta <- ta[rowSums(ta0[,-1:-3])!=0, ]
ta$Site <- factor(ta$Site, levels=c("Shitiping", "Jihuei", "Jialulan"), labels = c("A", "B", "C"))
ta$Season <- factor(ta$Season, levels = c("SP", "SU"), labels = c("4", "8"))

names(ta)[-1:-5] <- c("P. pangcahi", "Cyclopoapseudes sp.", "P. tagopilosus", "S. hansmuelleri", "I. multituberculata", "C. taitungensis", "P. setosa", "A. lenoprimorum","A. pedecerritulus", "T. nuwalianensis", "Z. shitipingensis", "Z. zorro")

ta$Microhabitat <- factor(ta$Microhabitat, 
                          levels = c("Amansia glomerata", "Asparagopsis sp.", "Asparagopsis taxiformis", "Caulerpa sp.", "Chlorodesmis sp.", "Corallina pilulifera", "Gelidium sp.", "Gravel, sand or silt", "Halimeda sp.", "Hypnea pannosa", "Hypnea sp.", "Jania sp.1", "Jania sp.2", "Mastophora rosea", "Padina sp.", "Plocamium sp.", "Tube of Eunice taoi"), 
                          labels = c("A. glomerata", "Asparagopsis sp.", "A. taxiformis", "Caulerpa sp.", "Chlorodesmis sp.", "C. pilulifera", "Gelidium sp.", "Gravel, sand, silt", "Halimeda sp.", "H. pannosa", "Hypnea sp.", "Jania sp.1", "Jania sp.2", "M. rosea", "Padina sp.", "Plocamium sp.", "Tube of E. taoi"))
ta$Microhabitat2 <- factor(ta$Microhabitat2, levels=c("Gravel, sand or silt", "Macroalgae", "Tube of Eunice taoi"), 
                           labels=c("Gravel, sand, silt", "Macroalgae", "Tube of E. taoi"))
names(ta)[2] <- "Month"

# Only keep data for analysis
b <- ta[, -1:-5]
rownames(b) <- with(ta, paste(Site, Month, Microhabitat, Replicate, sep="_"))

Preparing microhabitat grain size data

cs <- read_excel("../Data/Grain-size composition.xlsx", sheet = 6) %>% as.data.frame

cs$Microhabitat <- factor(cs$Microhabitat, levels=c("Amansia glomerata", "Asparagopsis sp.", "Asparagopsis taxiformis", 
                                                    "Caulerpa sp.", "Chlorodesmis sp.", "Corallina pilulifera", "Gelidium sp.",
                                                    "sand", "Halimeda sp.", "Hypnea pannosa", "Jania sp.1", "Jania sp.2", 
                                                    "Mastophora rosea", "Padina sp.",  "polychaete tube"), 
                          labels = c("A. glomerata", "Asparagopsis sp.", "A. taxiformis", "Caulerpa sp.", "Chlorodesmis sp.", 
                                     "C. pilulifera", "Gelidium sp.", "Gravel, sand, silt", "Halimeda sp.", "H. pannosa",  
                                     "Jania sp.1", "Jania sp.2", "M. rosea", "Padina sp.", "Tube of E. taoi"))
cs$Microhabitat2 <- as.character(cs$Microhabitat)
cs$Microhabitat2[!(cs$Microhabitat2=="Gravel, sand, silt" | cs$Microhabitat2=="Tube of E. taoi")] <- "Macroalgae"
cs$Microhabitat2 <- factor(cs$Microhabitat2)
cs$Site <- factor(cs$Site, levels=c("STP", "JH", "JLL"), labels=c("A", "B", "C"))
cs$Season <- factor(cs$Season, levels=c("SP", "SU"), labels=c("4", "8"))
cs <- cs[,c(1:3, 10, 4:9)]
names(cs)[2] <- "Month"

# Average by site, Month, and microhabitat
cs_avg <- aggregate(cs[,-1:-6], by=list(cs$Site, cs$Month, cs$Microhabitat, cs$Microhabitat2), FUN=mean, na.rm=TRUE)
names(cs_avg)[1:4] <- c("Site", "Month", "Microhabitat", "Microhabitat2")
cs_avg[is.na(cs_avg)] <- 0

cs[is.na(cs)] <- 0
cs_names <- with(cs,  paste(Site, Month, Microhabitat, Replicate, sep="_"))
ta_names <- with(ta,  paste(Site, Month, Microhabitat, Replicate, sep="_"))
env <- cs[match(ta_names, cs_names), -1:-6]
row.names(env) <- rownames(b)

# Logit transformation and normalized (divided by sd)
ens <- scale(logit(cs[, 7:10]/100))[match(ta_names, cs_names), ] %>% as.data.frame
row.names(ens) <- rownames(b)

Preparing CWB buoy data

cwb <-  fromJSON(file="../Data/C-B0050-001.json")

# Meta data 
meta <- lapply(cwb$cwbdata$Resources$Resource$Data$SeaSurfaceObs$Location, FUN=function(x)x[[1]])

# Monthly temperature
dat <- lapply(cwb$cwbdata$Resources$Resource$Data$SeaSurfaceObs$Location, FUN=function(x)x[[2]]$Monthly)

MonthlyTemperatureFun <- function(i){
  cbind(StationNameEN = lapply(meta, FUN=function(x)x$StationNameEN)[i] %>% as.character %>% as.factor,
        cbind(
          StationLatitude=lapply(meta, FUN=function(x)x$StationLatitude)[i] %>% as.numeric,
          StationLongitude=lapply(meta, FUN=function(x)x$StationLongitude)[i] %>% as.numeric,
          DataMonth=lapply(dat[[i]], FUN=function(x)x$DataMonth) %>% unlist %>% as.numeric,
          Maximum=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Maximum) %>% unlist %>% as.numeric,
          MaximumYear=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MaximumYear) %>% unlist %>% as.numeric,
          Mean=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Mean) %>% unlist %>% as.numeric,
          Minimum=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Minimum) %>% unlist %>% as.numeric,
          MinimumYear=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MinimumYear) %>% unlist %>% as.numeric,
          MinimumAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MinimumAnomaly) %>% unlist %>% as.numeric,
          MaximumAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MaximumAnomaly) %>% unlist %>% as.numeric,
          MeanAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MeanAnomaly) %>% unlist %>% as.numeric
          ) %>% as.data.frame
  )
}

# Hualien Buoy
hb <- MonthlyTemperatureFun(18) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")
# Chenggong
cg <- MonthlyTemperatureFun(8) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")
# Taitung Buoy
tb <- MonthlyTemperatureFun(34) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")

tmp <- rbind(hb, cg, tb)
tmp$Site <- factor(tmp$StationNameEN, levels=c("Hualien Buoy", "Chenggong", "Taitung Buoy"), labels = c("A", "B", "C"))
tmp$Month <- factor(tmp$DataMonth, levels=c(4, 8))

tmp_names <- with(tmp, paste(Site, Month, sep="_"))
ta_names <- with(ta,  paste(Site, Month, sep="_"))

env$Temperature <- tmp$Mean[match(ta_names, tmp_names)]
# Log 10 transformation and normalize (divided by sd)
ens$Temperature <- scale(log10(env$Temperature))

By replicates

Cluster Analysis

# Hierarchical Clustering based on square root transformed abundance data 
# Agglomeration used group average (= UPGMA) method
d <- vegdist(b^0.25)
hc <- hclust(d, method="average")
#plot(x = hc, labels =  row.names(hc), cex = 0.5, hang=-1)
#rect.hclust(tree = hc, k = 4)
fac <- ta[, 1:5]
fac$Group <- cutree(hc, 4)
# Reorder group from north to south as Groups 1, 2, 3
fac$Group <- factor(fac$Group, levels=c(1,2,3,4), labels = c(2, 3, 1, 4))

# Convert dendrogram to ggplot style
dhc <- as.dendrogram(hc)
ghc    <- dendro_data(dhc, type="rectangle") 

# Merge sample info to dendrogram label data frame
ord <- match(label(ghc)$label, rownames(b))
ghc[["labels"]]   <- cbind(ghc[["labels"]], fac[ord,])

# Group membership

# Group 1 (Shitiping = 82.6%)
table(subset(fac, Group==1)$Site)

## 
##  A  B  C 
## 19  2  2

# Group 2 (Jihuei = 52.6%, Jialulan = 28.9%)
table(subset(fac, Group==2)$Site)

## 
##  A  B  C 
##  7 20 11

# Group 3 (Jialulan = 77.6%, Jihuei = 17.2%)
table(subset(fac, Group==3)$Site)

## 
##  A  B  C 
##  3 10 45

Nonmetric Multidimensional Scaling

# Remove outlier JH_AP_Jania sp.2_1
md <- metaMDS(vegdist(b[-23,]^0.25))

## Run 0 stress 0.1395381 
## Run 1 stress 0.1381564 
## ... New best solution
## ... Procrustes: rmse 0.04435303  max resid 0.1245465 
## Run 2 stress 0.1406362 
## Run 3 stress 0.1401838 
## Run 4 stress 0.1484552 
## Run 5 stress 0.1420969 
## Run 6 stress 0.1474993 
## Run 7 stress 0.1407365 
## Run 8 stress 0.1477866 
## Run 9 stress 0.1380858 
## ... New best solution
## ... Procrustes: rmse 0.01438577  max resid 0.07559127 
## Run 10 stress 0.1399865 
## Run 11 stress 0.1456172 
## Run 12 stress 0.1415289 
## Run 13 stress 0.1463659 
## Run 14 stress 0.1380859 
## ... Procrustes: rmse 3.308223e-05  max resid 0.0001512511 
## ... Similar to previous best
## Run 15 stress 0.1398388 
## Run 16 stress 0.1419019 
## Run 17 stress 0.1434156 
## Run 18 stress 0.1405392 
## Run 19 stress 0.1404275 
## Run 20 stress 0.1423438 
## *** Best solution repeated 1 times

stress <- paste("Stress = ", deparse(round(md$stress,2)))

keep <- names(colSums(b[-23,])%>%sort(decreasing = TRUE))[1:5]

blab <- wascores(md$points, b[-23,]^0.25, expand = TRUE)[keep,] %>% as.data.frame
blab$label <- paste("italic(", sub(" ", "~~", rownames(blab)), ")", sep="")

# Only fits environmental vector to the replicate with grain size data
keep <- !env[-23, 1:4]%>%rowSums%>%is.na
set.seed(100)
(fit <- envfit(md$points[keep,], env[-23,][keep,]))

## 
## ***VECTORS
## 
##                       MDS1     MDS2     r2 Pr(>r)    
## Silt2Clay          0.29539  0.95538 0.0264  0.475    
## Fine2VeryFine      0.53423  0.84534 0.1250  0.031 *  
## Medium             1.00000  0.00269 0.2582  0.001 ***
## Coarse2VeryCoarse  0.88063 -0.47380 0.0912  0.078 .  
## Temperature       -0.69724  0.71684 0.0690  0.158    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

evec <- as.data.frame(scores(fit, display = "vectors")) * ordiArrowMul(fit, fill = 1.5)
evec$label <- row.names(evec)

# Mantel tests on each variables
man <- list()
man[[1]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Silt2Clay"]))
man[[2]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Fine2VeryFine"]))
man[[3]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Medium"]))
man[[4]] <- mantel(vegdist(b[-23,][keep,]^0.25), dist(env[-23,][keep, "Coarse2VeryCoarse"]))
man[[5]] <- mantel(vegdist(b[-23,]^0.25), dist(env[-23, "Temperature"]))

evec$mantel_r <- lapply(man, FUN=function(x)x$statistic) %>% unlist
evec$signif <- lapply(man, FUN=function(x)x$signif) %>% unlist
#evec$pvals <- cut(evec$signif, breaks=c(0, 0.01, 0.05, 0.1, 1), labels=c("<0.01", "<0.05", "<0.1", "<1")) %>% factor
evec$pvals <- factor(evec$signif)

print(evec)

##                         MDS1         MDS2             label    mantel_r signif
## Silt2Clay          0.1417114  0.458343873         Silt2Clay -0.05470974  0.827
## Fine2VeryFine      0.5576827  0.882441844     Fine2VeryFine  0.07029091  0.083
## Medium             1.5000000  0.004034451            Medium  0.21036458  0.001
## Coarse2VeryCoarse  0.7851920 -0.422451133 Coarse2VeryCoarse  0.11119698  0.006
## Temperature       -0.5406751  0.555876664       Temperature  0.13973555  0.001
##                   pvals
## Silt2Clay         0.827
## Fine2VeryFine     0.083
## Medium            0.001
## Coarse2VeryCoarse 0.006
## Temperature       0.001

# Merge MDS to environmental data frame
p2 <- ggplot(data=cbind(md$points, fac[-23,]), aes(x=MDS1, y=MDS2))+
  geom_point(alpha=1, size=5, stroke=1, aes(colour=Site, shape=Month))+
  geom_mark_hull(concavity = 10, aes(group=Group), colour="black", linetype=2)+
  geom_segment(data= evec[-5,], aes(x=0, y=0, xend=MDS1, yend=MDS2, linewidth=pvals, linetype=pvals),
               arrow = arrow(length=unit(.4, 'cm')), colour="red")+
  annotate("text", x=-0.9, y=1.1, label=stress, size=5) +
  annotate("text", x=1, y=0, label="1", size=10) +
  annotate("text", x=-0.5, y=-0.6, label="2", size=10) +
  annotate("text", x=-0.2, y=0.5, label="3", size=10) +
  annotate("text", x=1.7, y=1.2, label="b", size=7) +
  scale_shape_manual(values=c(1, 19))+
  scale_color_viridis_d()+
  scale_linewidth_manual(values=c(2, 1, 0.5, 0.5))+
  scale_linetype_manual(values=c(1, 1, 1, 2))+
  geom_label_repel(data=blab, aes(x=MDS1, y=MDS2, label=label), colour="black", 
             fill=gray(1, 0.6), size=4, label.padding = unit(0.25, "lines"), parse=T)+
  geom_label_repel(data=evec[-5,], aes(x=MDS1*1.2, y=MDS2*1.2, label=label), colour="red", 
             fill=gray(1, 0.6), size=4, label.padding = unit(0.25, "lines"), force=2)+
  labs(linewidth="p-value", linetype="p-value")+
  theme_bw() %+replace% large

ggdendrogram(hc, rotate = FALSE)+
  geom_point(data=label(ghc), aes(x, y, colour=Site, shape=Month), size=2, stroke=0.5)+
  geom_hline(yintercept = 0.8, linetype=2, colour="red")+
  annotate("text", x = 95, y = 0.85, label = "2", size = 5)+
  annotate("text", x = 42, y = 0.85, label = "3", size = 5)+
  annotate("text", x = 10, y = 0.85, label = "1", size = 5)+
  scale_shape_manual(values=c(1, 19))+
  scale_color_viridis_d()+
  labs(x= "Replicate",y = "Bray-Curtis Dissimilarity")+
  theme_bw() %+replace% theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), axis.line.y = element_line(colour="black")) %+replace%  rotate

m <- acast(dat=ta, Microhabitat~Site+Month, fun.aggregate = length) 
hc <- vegdist (m) %>% hclust(method="average")
clust_names <- row.names(m)[rev(hc$order)]
out <- cbind(md$points, fac[-23,])
out$Microhabitat <- factor(out$Microhabitat, levels=clust_names)

# Color code the microhabitat
cols <- c("#222222", "#f3c300", "#875692", "#f38400", "#a1caf1", "#be0032", rep("#c2b280", 7), rep("#848482", 3),"#008856")
shapes <- c(rep(19, 6), 0:6, 15, 17, 18:19)


# MDS for each site and Month
md <- lapply(splitBy(~Site+Month, ta), FUN=function(x)metaMDS(x[,-1:-5]))

## Wisconsin double standardization
## Run 0 stress 0.09177207 
## Run 1 stress 0.09177187 
## ... New best solution
## ... Procrustes: rmse 0.000338808  max resid 0.001136196 
## ... Similar to previous best
## Run 2 stress 0.09240116 
## Run 3 stress 0.1102589 
## Run 4 stress 0.09185124 
## ... Procrustes: rmse 0.07012022  max resid 0.1591274 
## Run 5 stress 0.09901149 
## Run 6 stress 0.09185037 
## ... Procrustes: rmse 0.06994872  max resid 0.1586776 
## Run 7 stress 0.1204544 
## Run 8 stress 0.09857996 
## Run 9 stress 0.09058073 
## ... New best solution
## ... Procrustes: rmse 0.1173281  max resid 0.306095 
## Run 10 stress 0.09240098 
## Run 11 stress 0.1198589 
## Run 12 stress 0.09185052 
## Run 13 stress 0.0955431 
## Run 14 stress 0.09050981 
## ... New best solution
## ... Procrustes: rmse 0.08364446  max resid 0.2192915 
## Run 15 stress 0.0909706 
## ... Procrustes: rmse 0.0747375  max resid 0.2755401 
## Run 16 stress 0.1253892 
## Run 17 stress 0.125494 
## Run 18 stress 0.09050979 
## ... New best solution
## ... Procrustes: rmse 0.0006690435  max resid 0.001675329 
## ... Similar to previous best
## Run 19 stress 0.09058075 
## ... Procrustes: rmse 0.08327585  max resid 0.2200089 
## Run 20 stress 0.09058272 
## ... Procrustes: rmse 0.08290963  max resid 0.2208044 
## *** Best solution repeated 1 times
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1345332 
## Run 1 stress 0.1480562 
## Run 2 stress 0.1468133 
## Run 3 stress 0.1551481 
## Run 4 stress 0.1383558 
## Run 5 stress 0.153372 
## Run 6 stress 0.1322342 
## ... New best solution
## ... Procrustes: rmse 0.09070836  max resid 0.298377 
## Run 7 stress 0.1403992 
## Run 8 stress 0.1755181 
## Run 9 stress 0.1678867 
## Run 10 stress 0.1564364 
## Run 11 stress 0.1442084 
## Run 12 stress 0.1524216 
## Run 13 stress 0.1361397 
## Run 14 stress 0.1403993 
## Run 15 stress 0.17642 
## Run 16 stress 0.153372 
## Run 17 stress 0.1346216 
## Run 18 stress 0.1678868 
## Run 19 stress 0.1550996 
## Run 20 stress 0.139706 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1208441 
## Run 1 stress 0.09427342 
## ... New best solution
## ... Procrustes: rmse 0.1234327  max resid 0.2720241 
## Run 2 stress 0.1098328 
## Run 3 stress 0.09427323 
## ... New best solution
## ... Procrustes: rmse 7.529209e-05  max resid 0.0002953634 
## ... Similar to previous best
## Run 4 stress 0.1382285 
## Run 5 stress 0.1349709 
## Run 6 stress 0.09332658 
## ... New best solution
## ... Procrustes: rmse 0.03245163  max resid 0.1489107 
## Run 7 stress 0.1196955 
## Run 8 stress 0.09427334 
## Run 9 stress 0.1212999 
## Run 10 stress 0.125197 
## Run 11 stress 0.1175232 
## Run 12 stress 0.1265351 
## Run 13 stress 0.1170811 
## Run 14 stress 0.09332658 
## ... New best solution
## ... Procrustes: rmse 0.000273151  max resid 0.0009519536 
## ... Similar to previous best
## Run 15 stress 0.1099111 
## Run 16 stress 0.09427323 
## Run 17 stress 0.09332666 
## ... Procrustes: rmse 0.0003160555  max resid 0.001152594 
## ... Similar to previous best
## Run 18 stress 0.1145061 
## Run 19 stress 0.1054994 
## Run 20 stress 0.1294716 
## *** Best solution repeated 2 times
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.138881 
## Run 1 stress 0.1380855 
## ... New best solution
## ... Procrustes: rmse 0.02290013  max resid 0.09354521 
## Run 2 stress 0.2044978 
## Run 3 stress 0.2363496 
## Run 4 stress 0.1750268 
## Run 5 stress 0.2121027 
## Run 6 stress 0.2069126 
## Run 7 stress 0.1518363 
## Run 8 stress 0.2125752 
## Run 9 stress 0.1393694 
## Run 10 stress 0.216466 
## Run 11 stress 0.1553091 
## Run 12 stress 0.133009 
## ... New best solution
## ... Procrustes: rmse 0.08724517  max resid 0.3626722 
## Run 13 stress 0.1638621 
## Run 14 stress 0.1400827 
## Run 15 stress 0.2068781 
## Run 16 stress 0.2026245 
## Run 17 stress 0.1564282 
## Run 18 stress 0.1526192 
## Run 19 stress 0.2163723 
## Run 20 stress 0.2174266 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     15: stress ratio > sratmax
##      5: scale factor of the gradient < sfgrmin
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.05720906 
## Run 1 stress 0.0572091 
## ... Procrustes: rmse 0.0002015069  max resid 0.0003602265 
## ... Similar to previous best
## Run 2 stress 0.0572092 
## ... Procrustes: rmse 0.0002634092  max resid 0.0004780251 
## ... Similar to previous best
## Run 3 stress 0.05720917 
## ... Procrustes: rmse 0.0002193702  max resid 0.0004674462 
## ... Similar to previous best
## Run 4 stress 0.05720909 
## ... Procrustes: rmse 7.852231e-05  max resid 0.00014762 
## ... Similar to previous best
## Run 5 stress 0.05720917 
## ... Procrustes: rmse 0.0002701912  max resid 0.0004919481 
## ... Similar to previous best
## Run 6 stress 0.05720914 
## ... Procrustes: rmse 0.0002475757  max resid 0.0004502463 
## ... Similar to previous best
## Run 7 stress 0.05720912 
## ... Procrustes: rmse 0.0002307228  max resid 0.0004188804 
## ... Similar to previous best
## Run 8 stress 0.08766255 
## Run 9 stress 0.08093205 
## Run 10 stress 0.08093206 
## Run 11 stress 0.1086705 
## Run 12 stress 0.05720904 
## ... New best solution
## ... Procrustes: rmse 4.555359e-05  max resid 7.944265e-05 
## ... Similar to previous best
## Run 13 stress 0.05720907 
## ... Procrustes: rmse 0.0001330564  max resid 0.0002378443 
## ... Similar to previous best
## Run 14 stress 0.08093209 
## Run 15 stress 0.05720903 
## ... New best solution
## ... Procrustes: rmse 3.619008e-05  max resid 0.0001015064 
## ... Similar to previous best
## Run 16 stress 0.05720904 
## ... Procrustes: rmse 2.007738e-05  max resid 5.202092e-05 
## ... Similar to previous best
## Run 17 stress 0.07810263 
## Run 18 stress 0.08093205 
## Run 19 stress 0.08093206 
## Run 20 stress 0.05720908 
## ... Procrustes: rmse 0.0001108149  max resid 0.0002103537 
## ... Similar to previous best
## *** Best solution repeated 3 times
## Wisconsin double standardization
## Run 0 stress 6.522281e-33 
## Run 1 stress 0 
## ... New best solution
## ... Procrustes: rmse 0.1307195  max resid 0.2126991 
## Run 2 stress 0 
## ... Procrustes: rmse 0.1368719  max resid 0.2230256 
## Run 3 stress 0 
## ... Procrustes: rmse 0.1450721  max resid 0.287975 
## Run 4 stress 0 
## ... Procrustes: rmse 0.1626683  max resid 0.2638379 
## Run 5 stress 0 
## ... Procrustes: rmse 0.09148311  max resid 0.1643084 
## Run 6 stress 1.910875e-05 
## ... Procrustes: rmse 0.149761  max resid 0.2413231 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1231937  max resid 0.2250026 
## Run 8 stress 0 
## ... Procrustes: rmse 0.1865359  max resid 0.3227937 
## Run 9 stress 0 
## ... Procrustes: rmse 0.1331486  max resid 0.2704025 
## Run 10 stress 0 
## ... Procrustes: rmse 0.115681  max resid 0.1553413 
## Run 11 stress 0 
## ... Procrustes: rmse 0.1302713  max resid 0.2177031 
## Run 12 stress 6.264849e-05 
## ... Procrustes: rmse 0.1032593  max resid 0.1609292 
## Run 13 stress 0 
## ... Procrustes: rmse 0.139783  max resid 0.2708793 
## Run 14 stress 0 
## ... Procrustes: rmse 0.1839108  max resid 0.3134858 
## Run 15 stress 0 
## ... Procrustes: rmse 0.1270628  max resid 0.1964893 
## Run 16 stress 0 
## ... Procrustes: rmse 0.172901  max resid 0.2816818 
## Run 17 stress 0 
## ... Procrustes: rmse 0.1588443  max resid 0.2648712 
## Run 18 stress 0 
## ... Procrustes: rmse 0.1501897  max resid 0.2424697 
## Run 19 stress 0 
## ... Procrustes: rmse 0.1613569  max resid 0.2847132 
## Run 20 stress 0 
## ... Procrustes: rmse 0.1396997  max resid 0.2310131 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin

# Stress levels
st <- lapply(md, FUN=function(x)x$stress) %>% ldply
st <- cbind(strsplit(st$.id, split="[|]") %>% ldply, st[, -1])
names(st) <- c("Site", "Month", "stress")
st$stress <- paste("Stress = ", round(st$stress,2))
st$Site <- factor(st$Site, levels=c("A", "B", "C"))
# MDS scores
md <- cbind(ta[,1:5], lapply(md, FUN=function(x)x$points) %>% ldply)
md$Microhabitat <- factor(md$Microhabitat, levels = clust_names)

ggplot(data=md, aes(x=MDS1, y=MDS2, colour=Microhabitat, shape=Microhabitat))+
  geom_point(size=3)+
  geom_text(data = st, inherit.aes = FALSE, aes(x=Inf, y = Inf, label = stress), vjust=1.2, hjust=1.05)+
  scale_colour_manual(values=cols)+
  scale_shape_manual(values=shapes)+
  facet_grid(Site~Month, scales="free")+
  theme_bw() %+replace% large

ggplot(data=md, aes(x=MDS1, y=MDS2, colour=Microhabitat2, shape=Microhabitat2))+
  geom_point(size=3)+
  geom_text(data = st, inherit.aes = FALSE, aes(x=Inf, y = Inf, label = stress), vjust=1.2, hjust=1.05)+
  scale_colour_manual(values=c("#222222", "#f3c300", "#875692"))+
  scale_shape_manual(values=c(19, 2, 15))+
  labs(colour="Microhabitat", shape="Microhabitat")+
  facet_grid(Site~Month, scales="free")+
  theme_bw() %+replace% large

drows <- vegdist(t(b)^0.25)
dcols <- vegdist(b^0.25)

annotation_col <- fac[, -4:-5][,c(3,2,1,4)]
rownames(annotation_col) <- row.names(b)

annotation_colors <- list(
   Microhabitat=c("Tube of E. taoi" = "#222222", "A. glomerata" = "#f3c300", 
                  "M. rosea" = "#875692", "Jania sp.1" = "#f38400", 
                  "Gravel, sand, silt" = "#a1caf1","H. pannosa" = "#be0032", 
                  "A. taxiformis" = "#c2b280", "Asparagopsis sp." = "#c2b280", 
                  "Chlorodesmis sp." = "#c2b280", "Halimeda sp." = "#c2b280", 
                  "Padina sp." = "#c2b280", "C. pilulifera" = "#c2b280", 
                  "Caulerpa sp." = "#c2b280", "Plocamium sp." = "#848482", 
                  "Hypnea sp." = "#848482", "Gelidium sp." = "#848482", 
                  "Jania sp.2" = "#008856"),
  Month=c("4"="white", "8"="black"),
  Site=c(A="#440154FF", B="#21908CFF", C="#FDE725FF"),
  Group=c("1"="#440154FF", "2"="#21908CFF", "3"="#FDE725FF", "4"="white")
)

pheatmap(t(b)^0.25, 
         color = brewer.pal(9, "Blues"), 
         clustering_distance_rows = drows, 
         clustering_distance_cols = dcols, 
         clustering_method = "average", 
         legend_breaks = 0:4,
         legend_labels = (0:4)^4,
         cutree_cols = 4, 
         annotation_col = annotation_col, 
         annotation_colors = annotation_colors,
         fontsize_col=4, fontsize_row=12)

annotation_col <- fac[, c(-3, -5)][,c(3,2,1,4)]
rownames(annotation_col) <- row.names(b)
names(annotation_col)[1] <- "Microhabitat"

annotation_colors <- list(
  Microhabitat=c("Tube of E. taoi" = "#222222", "Gravel, sand, silt" = "#a1caf1", "Macroalgae" = "#c2b280"),
  Month=c("4"="white", "8"="black"),
  Site=c(A="#440154FF", B="#21908CFF", C="#FDE725FF"),
  Group=c("1"="#440154FF", "2"="#21908CFF", "3"="#FDE725FF", "4"="white")
)

pheatmap(t(b)^0.25, 
         color = brewer.pal(9, "Blues"), 
         clustering_distance_rows = drows, 
         clustering_distance_cols = dcols, 
         clustering_method = "average", 
         legend_breaks = 0:4,
         legend_labels = (0:4)^4,
         cutree_cols = 4, 
         annotation_col = annotation_col, 
         annotation_colors = annotation_colors,
         fontsize_col=4, fontsize_row=12)

Aggregate by site-Month

ta3 <- aggregate(b, by = list(ta$Site, ta$Month), FUN=mean)
names(ta3)[1:2] <- c("Site", "Month")
b3 <- ta3[, -1:-2]
rownames(b3) <- with(ta3, paste(Site, Month, sep="_"))

# Average by site and Month
cs_avg2 <- aggregate(cs[,-1:-6], by=list(cs$Site, cs$Month), FUN=mean, na.rm=TRUE)
names(cs_avg2)[1:2] <- c("Site", "Month")
cs_avg2

##   Site Month  Silt2Clay Fine2VeryFine    Medium Coarse2VeryCoarse
## 1    A     4 0.07468493      9.590777 12.213665          42.40659
## 2    B     4 0.54811731      8.016581 10.550036          55.88527
## 3    C     4 0.72627053     12.150971  3.015219          12.67897
## 4    A     8 2.85256410     16.025641 18.878205          62.24359
## 5    B     8 0.93316520     10.156021  9.265062          59.64575
## 6    C     8 0.65368487     15.726078  5.927816          33.94242

Cluster Analysis

# Hierarchical Clustering based on square root transformed abundance data 
# Agglomeration used group average (= UPGMA) method
(d <- vegdist(b3^0.25))

##           A_4       B_4       C_4       A_8       B_8
## B_4 0.3183515                                        
## C_4 0.3953006 0.3404144                              
## A_8 0.3071157 0.3873459 0.4548924                    
## B_8 0.3853953 0.1707948 0.2189590 0.4929781          
## C_8 0.3759794 0.1824998 0.2880368 0.3637466 0.1863058

hc <- hclust(d, method="average")
fac <- ta3[, 1:2]

# Convert dendrogram to ggplot style
dhc <- as.dendrogram(hc)
ghc    <- dendro_data(dhc, type="rectangle") 

# Merge sample info to dendrogram label data frame
ord <- match(label(ghc)$label, rownames(b3))
ghc[["labels"]]   <- cbind(ghc[["labels"]], fac[ord,])

p5 <- ggdendrogram(hc, rotate = FALSE)+
  geom_point(data=label(ghc), aes(x, y, colour=Site, shape=Month), size=4, stroke=1.3)+
  annotate("text", x=6, y=0.4, label="a", size=7) +
  scale_shape_manual(values=c(1, 19))+
  scale_color_viridis_d()+
  labs(x= "",y = "Bray-Curtis Dissimilarity")+
  theme_bw() %+replace% theme(legend.position = "none", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(),axis.line.y = element_line(colour="black")) %+replace% large

Comibining site-season clustering with replicate-level nMDS

p5/p2+plot_layout(nrow = 2, heights = c(1, 3))